/*
 * DiagonalPreconditioner.cs
 * 
 * Copyright (c) 2003-2004, dnAnalytics Project. All rights reserved.
*/

using System;
using dnAnalytics.Exceptions;
using dnAnalytics.Resources;

namespace dnAnalytics.LinearAlgebra.Solvers.Preconditioners
{
    /// <summary>
    /// A diagonal preconditioner. The preconditioner uses the inverse
    /// of the matrix diagonal as preconditioning values.
    /// </summary>
    /// <include file='../../../../../examples/trunk/LinearAlgebra/Solvers/Preconditioners/DiagonalPreconditioner.xml' path='example'/> 
    public sealed class DiagonalPreconditioner : IPreconditioner
    {
        /// <summary>
        /// The coefficient matrix upon which this preconditioner is based.
        /// </summary>
        private SparseMatrix mCoefficientMatrix;
        /// <summary>
        /// The inverse of the matrix diagonal.
        /// </summary>
        private double[] mInverseDiagonals;
        /// <summary>
        /// A flag used to indicate if the preconditioner has been created, i.e.
        /// if the <b>mInverseDiagonals</b> array has been filled with the correct
        /// values.
        /// </summary>
        private bool mWasPreconditionerCreated;

        /// <summary>
        /// Creates a <b>DiagonalPreconditioner</b> instance.
        /// </summary>
        /// <param name="coefficientMatrix">The <see cref="SparseMatrix"/> upon which this
        /// preconditioner is based.</param>
        /// <exception cref="ArgumentNullException">If <paramref name="coefficientMatrix"/> is <b>null</b>.</exception>
        /// <exception cref="NotSquareMatrixException">If <paramref name="coefficientMatrix"/> is not a square matrix.</exception>
        public DiagonalPreconditioner(Matrix coefficientMatrix)
        {
            if (coefficientMatrix == null)
            {
                throw new ArgumentNullException("coefficientMatrix", Strings.NullParameterException);
            }

            if (coefficientMatrix.Rows != coefficientMatrix.Columns)
            {
                throw new NotSquareMatrixException();
            }

            mCoefficientMatrix = new SparseMatrix(coefficientMatrix);
        }

        /// <summary>
        /// Initializes the preconditioner and loads the internal data structures.
        /// </summary>
        /// <remarks>
        /// The preconditioner is created by copying the diagonals of the original 
        /// matrix and then storing the inverses of these diagonal values.
        /// </remarks>
        public void CreatePreconditioner()
        {
            if (!mWasPreconditionerCreated)
            {
                mInverseDiagonals = new double[mCoefficientMatrix.Rows];
                for (int i = 0; i < mCoefficientMatrix.Rows; i++)
                {
                    mInverseDiagonals[i] = 1 / mCoefficientMatrix[i, i];
                }
                mWasPreconditionerCreated = true;
            }
        }

        /// <summary>
        /// Solves the matrix equation <b>Ax = b</b>.
        /// </summary>
        /// <param name="vector">The right hand side vector.</param>
        /// <param name="result">The left hand side vector.</param>
        /// <exception cref="ArgumentNullException">If <paramref name="vector"/> is <b>null</b>.</exception>
        /// <exception cref="ArgumentNullException">If <paramref name="result"/> is <b>null</b>.</exception>
        /// <exception cref="NotConformableException"><para>If <paramref name="vector"/> and <paramref name="result"/> do not
        /// have the same size.</para>
        /// <para> - or -</para>
        /// <para>If the size of <paramref name="result"/> is different the number of rows of the coefficient matrix.</para></exception>
        public void Solve(Vector vector, Vector result)
        {
            // Solve Dx = y
            // Since we stored the inverses of the diagonal values we can do:
            // x = D^-1 * y

            if (vector == null)
            {
                throw new ArgumentNullException("vector", Strings.NullParameterException);
            }

            if (result == null)
            {
                throw new ArgumentNullException("result", Strings.NullParameterException);
            }

            if ((result.Count != vector.Count) || (result.Count != mCoefficientMatrix.Rows))
            {
                throw new NotConformableException(Strings.ParametersNotConformable);
            }

            if (!mWasPreconditionerCreated)
            {
                CreatePreconditioner();
            }

            for (int i = 0; i < mInverseDiagonals.Length; i++)
            {
                result[i] = vector[i] * mInverseDiagonals[i];
            }
        }

        /// <summary>
        /// Solves the matrix equation <b>Ax = b</b>.
        /// </summary>
        /// <remarks>
        /// When using the <b>Solve</b> method consider carefully if the
        /// left hand side vector can be reused. If so it would be better to 
        /// use the overloaded <b>Solve</b> method.
        /// </remarks>
        /// <param name="vector">The right hand side vector.</param>
        /// <returns>The left hand side vector.</returns>
        /// <exception cref="ArgumentNullException">If <paramref name="vector"/> is <b>null</b>.</exception>
        /// <exception cref="NotConformableException"><para>If the size of <paramref name="vector"/> is different the number of rows of the coefficient matrix.</para></exception>
        public Vector Solve(Vector vector)
        {
            if (vector == null)
            {
                throw new ArgumentNullException("vector", Strings.NullParameterException);
            }

            if (vector.Count != mCoefficientMatrix.Columns)
            {
                throw new NotConformableException(Strings.ParameterNotConformable);
            }

            Vector result = VectorBuilder.CreateVector(vector.Count, VectorType.Dense);
            Solve(vector, result);
            return result;
        }
    }
}
